A Practical Primer in Human Complex Genetics

with a use-case in cardiovascular disease

1 About this primer

Ever since the first genome-wide association study (GWAS) on age-related macular degeneration, and the promise of personalized medicine in the wake of the Human Genome Project, large-scale genetic association studies hold significant sway in contemporary health research and drive drug-development pipelines. In the past 2 decades, researchers delved into GWAS, aiming to unveil genetic variations linked to both human traits, such as the color of your eyes, and rare and common complex diseases. These findings serve as crucial keys to unravel the intricate mechanisms underlying diseases, shedding light on whether the correlations identified in observational studies between risk factors and diseases are truly causal.

These studies have ushered in an exciting era where many researchers thrive on developing new methods and bioinformatic tools to parse ever-growing large datasets collected large population-based biobanks. However, the analyses of these data are challenging and it can be daunting to see the forest for tree among the many tools and their various functions. Enter A Practical Primer in Human Complex Genetics. This GitBook was originally written back in 2022 for the Genetic Epidemiology course organized by the Master Epidemiology of Utrecht University. This practical guide will teach you how to design a GWAS, perform quality control (QC), execute the actual analyses, annotate the GWAS results, and perform further downstream post-GWAS analyses. Throughout the book you’ll work with ‘dummy’, that is fake, data, but in the end, we will use real-world data from the first release of the Welcome Trust Case-Control Consortium (WTCCC) focusing on coronary artery disease (CAD).

A major component of modern-day GWAS is genetic imputation, but for practical reasons it is not part of this book. However, I will provide some pointers as to how to go about do this with minimal coding or scripting experience. Likewise, the courses does not cover the aspects of meta-analyses of GWAS, but some excellent resources exist to which I will direct. As this practical primer evolves, these and other topics may find their place in this book. I should also point out that emphasis of this book is on it being a practical primer. It is intended to provide some practical guidance to doing GWAS, and while theory is important, I will not cover this. Again, some very useful and excellent work exists to which I will point you, but I really want you to learn - and understand the theory - by doing.

So, although originally crafted as a companion for the course, this practical guide stands on its own as a comprehensive resource for diving into all facets of doing a GWAS — save for experimental follow-up, of course 😉.

I can imagine this seems overwhelming, but trust me, you’ll be okay. Just follow this practical. You’ll learn by doing and at the end of the day, you can execute a GWAS independently.

Ready to start?

2 Some background reading

Standing on the shoulders of giants, that’s what this book and I do. I want to acknowledge some great work that has helped me tremendously and, really, this book wouldn’t exist without this awesome work. So, I do want to give you some background reading. Is it a prerequisite? No, not really. For starters, the course covers most and you’ll learn as you go. And if you didn’t come here through the course, you’ll be fine just the same. That said, it’s a always good idea to get familiar with these works as you move forward on your path towards your first GWAS - in fact, I had these printed out with markings and writings all over them as I executed my first GWAS, and they’ve been great as a reference many times after.

Large parts of this work are based on four awesome Nature Protocols from the Zondervan group at the Wellcome Center Human Genetics.

1. Zondervan KT et al. Designing candidate gene and genome-wide case-control association studies. Nat Protoc 2007.
2. Pettersson FH et al. Marker selection for genetic case-control association studies. Nat Protoc 2009.
3. Anderson CA et al. Data QC in genetic case-control association studies. Nat Protoc 2010.
4. Clarke GM et al. Basic statistical analysis in genetic case-control studies. Nat Protoc 2011.

An update on the community standards of QC for GWAS can be found here:

1. Laurie CC et al. Quality control and quality assurance in genotypic data for genome-wide association studies. Genet Epidemiol 2010.

With respect to imputation and meta-analyses of GWAS you should also get familiar with the following two works:

1. Marchini, J. and Howie, B. Genotype imputation for genome-wide association studies. Nat Rev Genet 2010
2. de Bakker PIW et al. Practical aspects of imputation-driven meta-analysis of genome-wide association studies. Hum Mol Genet 2008.
3. Winkler TW et al. Quality control and conduct of genome-wide association meta-analyses. Nat Protoc 2014.

Are you ready?

Are you ready? Did you bring coffee and a good dose of energy? Let’s start! Your first point of action is to prepare your system for this course in Chapter 3.

3 Getting started

Before getting started, we need to discuss computers, specifically your computer. Most programs made to execute genetic epidemiology studies are developed for the Unix environment, for example Linux and macOS. So, they may not work as intended in a Windows environment. Windows does allow users to install a linux subsystem within Windows 10+ and you can find the detail guide here.

However, I highly recommend one of two options.

• One, install a linux subsystem on your Windows computer (for example a virtual machine with Ubuntu could work).
• Two, switch to macOS in combination with homebrew. This will give you all the flexibility to use Unix-based programs for your genetic epidemiology work and at the same time you’ll keep the advantage of a powerful computer with a user-friendly interface.

I chose the latter.

For this practical every command is intended for Linux/macOS, in other words Unix-systems.

3.1 CoCalc vs. Standalone

For the purpose of this practical primer there are one of two steps you can to take to get started.

• When you are following the course, you were invited to the CoCalc server of this course. So, you will want to read the section CoCalc.
• When you want to use this book as a standalone, you should check out the instructions in section Standalone - this is probably also the section you want to follow for real-world cases.

But first, I’ll briefly provide some background on the various programs that are commonly used Chapter 4.

4 Getting started - programs

We’ll use a few programs throughout this practical. You’ll probably need these for your (future) genetic epidemiology work too (Table 4.1).

Table 4.1: Programs needed for genetic epidemiology.

Program	Link	Description
PLINK	https://www.cog-genomics.org/plink2/	PLINK is a free, open-source genetic analysis tool set, designed to perform a range of basic data parsing and quality control, as well as basic and large-scale analyses in a computationally efficient manner.
R	https://cran.r-project.org/	A program to perform statistical analysis and visualizations.
RStudio	https://www.rstudio.com	A user-friendly R-wrap-around for code editing, debugging, analyses, and visualization.
Homebrew	https://brew.sh	A great extension for Mac-users to install really useful programs that Apple didn't.

4.1 RStudio

RStudio is a very user-friendly interface around R that makes your R-scripting-life a lot easier. You should get used to that. RStudio comes with R so you don’t have to worry about that.

4.2 PLINK

Right, onto PLINK.

All genetic analyses can be done in PLINK, even on your laptop, but with large datasets, for example UK Biobank size, it is better to switch to a high-performance computing cluster (HPC) like we have available at the Utrecht Science Park. The original PLINK v1.07 can be found here, but nowadays we are using a newer, faster version: PLINK v1.9 which can be found here. It still says ‘PLINK 1.90 beta’ (Figure 4.1), but you can consider this version stable and save to work with, but as you can see, some functions are not supported anymore.

Figure 4.1: The PLINK v1.9 website.

4.2.1 Alternatives to PLINK

Nowadays, a lot of people also use programs like SNPTEST, BOLT-LMM, GCTA, or regenie as alternatives to execute GWAS. These programs were designed with specific use-cases in mind, for instance really large biobank data including hundreds of thousands individuals, better control for population stratification, the ability to estimate trait heritability or Fst, and so on.

4.3 Other programs

Mendelian randomization can be done either with the SMR or GSMR function from GCTA, or with R-packages, like TwoSampleMR.

4.4 Did you choose?

Now is the moment of decision. You should opt for the CoCalc-instructions when you following the course. In all other cases, your best bet are the Standalone-instructions.

Go to Chapter 5 for the Standalone-instructions.

Go to Chapter 6 for the CoCalc-instructions.

5 Getting started - Standalone

So, you plan to use this book as ‘Standalone’ on a macOS environment. This means you’ll need to install a few things first.

5.1 The data you need

You’ll need to start by downloading the data you need for this practical to your Desktop.

Here’s the link to the data.

Link to Google Drive with data

Make sure you put the data in the ~/Desktop/practical/ folder.

The data are pretty large (approx. 15Gb), so this will take a minute or two depending on your internet connection. Time to stretch your legs or grab a coffee (data scientists don’t drink tea).

Note: the WTCCC1 data are not available for public download, as you need permission to use them. However, we do provide the instructions under the assumption your WTCCC1 data will be on genome build 37 and imputed against 1000G phase 3.

5.2 Terminal

For all the programs we use, except RStudio, you will need the Terminal. This comes with every major operating system; on Windows it is called ‘PowerShell’, but let’s not go there. And regardless, you will (have to start to) make your own scripts. The benefit of using scripts is that each step in your workflow is clearly stipulated and annotated, and it allows for greater reproducibility, easier troubleshooting, and scaling up to high-performance computer clusters.

Open the Terminal, it should be on the left in the toolbar as a little black computer-monitor-like icon. Mac users can type command + space and type terminal, a Terminal screen should open.

From now on we will use little code blocks like the example to indicate a code you should type/copy-paste and hit enter. If a code is followed by a comment, it is indicated by a # - you don’t need to copy-paste and execute this.

CODE BLOCK

CODE BLOCK # some comment here

5.3 Navigating the Terminal

You can navigate around the computer through the terminal by typing cd <path>; cd stands for “change directory” and <path> means “some_file_directory_you_want_to_go_to”.

This command will bring you to your home directory.

cd ~ 

This will bring you to the parent directory (up one level).

cd ../ 

This will bring you to the XXX directory.

cd XXX 

Let’s navigate to the folder you just downloaded.

cd ~/Desktop/practical

Let’s check out what is inside the directory, by listing (ls) its contents.

This command shows files as list; the -l makes it a vertical list and adds more information, you can also remove it and simply type ls - go on, and try.

ls -l 

This command shows files as list with human readable format.

ls -lh 

Adding the flags -lh will get you the contents of a directory in a list (-l) and make the size ‘human-readable’ (-h).

Adding -t shows the files as list sorted by time edited.

ls -lt 

Adding -S shows the files as list sorted by size.

ls -lS 

You can also count the number of files. Just ‘pipe’ the result from ls to the next program wc (‘wordcount’) and list the number of lines, -l. In this case -l is a flag used by wc and it has a different meaning than it does for ls.

ls | wc -l

And if you want to know all the function of a program simply type the following.

man ls

This will take you to a manual of the program with an extensive description of each flag (Figure 5.1).

Figure 5.1: Partial output from the ls manual.

5.4 Installing the software

5.4.1 brew

Linux has a great package-manager that is lacking on macOS. You can install brew to compensate for this. This adds the ability to install almost any Linux-based program through the Terminal such as wget, llvm, etc.

Open Terminal and execute the following:

/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

Check if everything is in order.

brew doctor

It shouldn’t report any errors.

5.4.2 PLINK

First, we’ll get PLINK. Navigate to the PLINK v1.9 website, which can be found here. Download the macOS (64-bit) version under ‘Stable (beta x.x, day month year)’.

Note: Apple produced Intel-based computers for a few years back, and most programs, packages, libraries and whatnot are designed for that. So, I highly recommend using software designed for that and activating Rosetta2 in your Terminal. Don’t know how to do that? Following these instructions.

Unzip the folder and put plink in the practical folder.

mv -v ~/Downloads/plink_mac_20231211/plink ~/Desktop/practical/plink 

5.4.3 Installing R and RStudio

Let’s go ahead and use brew to install the R and RStudio software.

In Terminal execute the following and just follow the instructions.

brew install rstudio
brew install --cask r

Now close the terminal window - really make sure that the terminal-program has quit.

Open your fresh installation of RStudio by double clicking the icon. You should be seeing something like figure 5.2

Figure 5.2: RStudio screenshot.

In the top right, you see a little green-white plus-sign, click this and select ‘R Notebook’ (Figure 5.3).

Figure 5.3: RStudio screenshot.

You will create an untitled (Untitled1) R notebook: you can combine text descriptions, like you would in a lab-journal, with code-sections. Read what is in the notebook to get a grasp on that (Figure 5.4).

Figure 5.4: RStudio screenshot.

Right, you should be installing some packages. To do so, you can remove plot(cars) (or leave and create a new code-block as per instructions in the notebook), and copy paste the code below. Make sure to put in a code block like the example in which plot(cars) is in.

remotes::install_github(c("rstudio/rmarkdown"))

install.packages(c("formatR", "remotes", 
                   "httr", "usethis", 
                   "data.table", "devtools", 
                   "dplyr", "tibble", "tidyverse", 
                   "openxlsx",
                   "ggplot2",
                   "ggsci", "ggthemes",
                   "qqman", "CMplot", "plotly", 
                   "openxlsx"))
devtools::install_github("kassambara/ggpubr")

devtools::install_github("oliviasabik/RACER")

remotes::install_github("MRCIEU/TwoSampleMR")
devtools::install_github("MRCIEU/MRInstruments")

if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("geneplotter")

You should load these packages too.

library(rmarkdown)
library(formatR)

library(openxlsx)

library(data.table)

library(tibble)
library(tidyverse)
library(dplyr)
library(plotly)

library(ggplot2)
library(devtools)
library(ggpubr)
library(ggsci)
library(ggthemes)

library(qqman)
library(CMplot)
library(RACER)

library(remotes)
library(TwoSampleMR)
library(MRInstruments)

library("geneplotter")

All in all this may take some time, good moment to relax, review your notes, stretch your legs, or take a coffee.

5.5 Are you ready?

Are you ready? Did you bring coffee and a good dose of energy? Let’s start!

Oh, one more thing: you can save your notebook, the one you just created, to keep all the R codes you are applying in the next chapters and add descriptions and notes. If you save this notebook you’ll notice that a html-file is created. This file is a legible webbrowser-friendly version of your work and contains the codes and the output (code messages, tables, and figures). And the nice thing is, that you can easily share it with others over email.

Ok. ’Nough said, let’s move on to cover some basics in Chapter 7.

6 Getting started - CoCalc

Now, pay attention. If you came here through the course Genetic Epidemiology, you don’t have to do anything. All the data you need are already downloaded and everything was set up for on the server, just follow the link as provided by the course-instructors.

Here we provide a few short instructions to navigate CoCalc.

6.1 Starting the course on CoCalc

Once logged in you should see a screen similar to the one below.

Figure 6.1: CoCalc after logging in.

Navigate to the Course by clicking on the blue Your UBC CoCalc 1 projects…-button and selecting the course Genetic Epidemiology you are following.

Figure 6.2: CoCalc projects.

Once you are in the course, you will see a screen like this one below.

Figure 6.3: CoCalc student page.

In the Course Book you’ll find the handouts you need. In the share data folder you will find the data you’ll need.

6.2 Beginning your own notebook

If you want to start your own notebook, you can do so by clicking on the New button and selecting Jupyter Notebook.

(#fig:cocalc_notebook)CoCalc new notebook.

You will probably get a screen asking you to select a kernel. You can choose R.

Figure 6.4: CoCalc kernel selection.

Now you are ready to start your own notebook - like below. You can create a new ‘cell’ with format ‘code’ or ‘text’ or ‘markdown’ and start typing.

Figure 6.5: CoCalc starting your notebook.

6.3 Are you ready?

Are you ready? Did you bring coffee and a good dose of energy? Let’s start!

Oh, one more thing: you can save your notebook, the one you just created, to keep all the R codes you are applying in the next chapters and add descriptions and notes. If you save this notebook you’ll notice that a html-file is created. This file is a legible webbrowser-friendly version of your work and contains the codes and the output (code messages, tables, and figures). And the nice thing is, that you can easily share it with others over email.

Ok. ’Nough said, let’s move on to cover some basics in Chapter 7.

7 Steps in a Genome-Wide Association Study

Now that you understand a bit of the navigation in Unix-systems, you’re ready to start this practical primer. We will make use of a dummy dataset containing cases and controls. We will explain and execute the following steps:

1. convert raw data to a more memory-efficient format
2. apply extensive quality control on samples and SNPs
3. assess the ancestral background of your study population
4. perform association testing
5. visualize association results

7.1 Converting datasets

The format in which genotype data are returned to investigators varies among genome-wide SNP platforms and genotyping centers. Usually genotypes have been called by a genotyping center and returned in the standard PED and MAP file formats designed for PLINK.

A PED file is a white space (space or tab)-delimited file in which each line represents one individual and the first six columns are mandatory and in the following order:

• ‘Family ID’,
• ‘Individual ID’,
• ‘Paternal ID’,
• ‘Maternal ID’,
• ‘Sex (1=male, 2=female, 0=missing)’, and
• ‘Phenotype (1=unaffected, 2=affected, 0=missing)’.

The subsequent columns denote genotypes that can be any character (e.g., 1, 2, 3, 4 or A, C, G, T). Zero denotes a missing genotype. Each SNP must have two alleles (i.e., both alleles are either present or absent). The order of SNPs in the PED file is given in the MAP file, in which each line denotes a single marker and the four white-space–separated columns are chromosome (1–22, X, Y or 0 for unplaced), marker name (typically an rs number), genetic distance in Morgans (this can be fixed to 0) and base-pair position (bp units).

Let’s start by using PLINK to converting the datasets to a lighter, binary form (a .bed-file). This file saves data in a more memory- and time-efficient manner (in a ‘binary’-format) to facilitate the analysis of large-scale data sets (Purcell S. 2007). The marker-information is stored in the .bim-file and the family information in the .fam-file. PLINK creates a .log file (named raw-GWA-data.log) that details (among other information) the implemented commands, the number of cases and controls in the input files, any excluded data and the genotyping rate in the remaining data. This file is very useful for checking whether the software is successfully completing commands.

Make sure you are in the right directory. Do you remember how to get there?

cd ~/Desktop/practical

Next, we’ll make a project directory.

mkdir -v ~/Desktop/practical/dummy_project

Now, we’ll convert the .ped/.map files to the binary-format.

plink --file rawdata/raw-GWA-data --make-bed --out dummy_project/rawdata

Let’s review the .log-file for a bit. It should look something like this:

PLINK v1.90b7.2 64-bit (11 Dec 2023)           www.cog-genomics.org/plink/1.9/
(C) 2005-2023 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to dummy_project/rawdata.log.
Options in effect:
  --file rawdata/raw-GWA-data
  --make-bed
  --out dummy_project/rawdata

16384 MB RAM detected; reserving 8192 MB for main workspace.
.ped scan complete (for binary autoconversion).
Performing single-pass .bed write (317503 variants, 2000 people).
--file: dummy_project/rawdata-temporary.bed +
dummy_project/rawdata-temporary.bim + dummy_project/rawdata-temporary.fam
written.
317503 variants loaded from .bim file.
2000 people (997 males, 1003 females) loaded from .fam.
2000 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 2000 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Warning: 11440 het. haploid genotypes present (see dummy_project/rawdata.hh );
many commands treat these as missing.
Total genotyping rate is 0.985682.
317503 variants and 2000 people pass filters and QC.
Among remaining phenotypes, 1023 are cases and 977 are controls.
--make-bed to dummy_project/rawdata.bed + dummy_project/rawdata.bim +
dummy_project/rawdata.fam ... done.

So, there are 317,503 variants included for 2,000 people, 997 males and 1,003 females. All of these individuals are ‘founders’. There are 1,023 cases and 977 controls. The genotyping rate is about 98.6% which is pretty good. Lastly, there are 11,440 heterozygous haploid genotypes present.

Question: Can you think off what the ‘11,440 heterozygous haploid genotypes present’ represent?

7.2 Quality control

We are ready for some quality control and quality assurance, heavily inspired by Anderson et al. (Anderson C. A. 2010) and Laurie et al. (Laurie C. C. 2010). In general, we should check out a couple of things regarding the data quality on two levels:

1. samples
2. variants

So, we will investigate the following:

• Are the sexes based on genetic data matching the ones given by the phenotype file?
• Identify individuals that are outliers in terms of missing data (call rate) or heterozygosity rates. This could indicate a genotyping error or sample swap.
• Identify duplicated or related individuals.
• Identify individuals with divergent ancestry.
• What are the allele frequencies?
• What is the per-SNP call rate?
• In the case of a case-control study (which is the case here), we need to check differential missingness between cases and controls.

Question: Can you think of other scenarios in which you may want to extend the check on differential missingness beyond a check between cases and controls?

7.3 Let’s get our hands dirty

All clear? Let’s start the work. On to step 1 of the QC for GWAS: filter samples of poor quality in Chapter 8.

8 Sample QC

Let’s start with the per-sample quality control. The goal here is to identify all the samples that for some reason (mentioned in the previous section) should be excluded.

8.1 Sex

We need to identify of individuals with discordant sex information comparing phenotypic and genotypic data. Let’s calculate the mean homozygosity rate across X-chromosome markers for each individual in the study.

plink --bfile dummy_project/rawdata --check-sex --out dummy_project/rawdata

This produces a file with the following columns:

• FID Family ID
• IID Within-family ID
• PEDSEX Sex code in input file
• SNPSEX Imputed sex code (1 = male, 2 = female, 0 = unknown)
• STATUS ‘OK’ if PEDSEX and SNPSEX match and are nonzero, ‘PROBLEM’ otherwise
• F Inbreeding coefficient, considering only X chromosome. Not present with ‘y-only’.
• YCOUNT Number of nonmissing genotype calls on Y chromosome. Requires ‘ycount’/‘y-only’.

We need to get a list of individuals with discordant sex data.

cat dummy_project/rawdata.sexcheck | awk '$5 =="STATUS" || $5 =="PROBLEM"'  > dummy_project/rawdata.sexprobs.txt

Let’s have a look at the results.

cat dummy_project/rawdata.sexprobs.txt

Table 8.1: Sex issues.

FID	IID	PEDSEX	SNPSEX	STATUS	F
772	772	2	0	PROBLEM	0.3084
853	853	2	0	PROBLEM	0.3666
1,920	1,920	2	0	PROBLEM	0.4066

When the homozygosity rate (F) is more than 0.2, but less than 0.8, the genotype data are inconclusive regarding the sex of an individual and these are marked in column SNPSEX with a 0, and the column STATUS “PROBLEM”.

Report the IDs of individuals with discordant sex information (Table 8.1) to those who conducted sex phenotyping. In situations in which discrepancy cannot be resolved, add the family ID (FID) and individual ID (IID) of the samples to a file named fail-sexcheck-qc.txt (one individual per line, tab delimited).

grep "PROBLEM" dummy_project/rawdata.sexcheck | awk '{ print $1, $2}'  > dummy_project/fail-sexcheck-qc.txt

8.2 Sample call rates

Let’s get an overview of the missing data per sample and per SNP.

plink --bfile dummy_project/rawdata --missing --out dummy_project/rawdata

This produces two files, rawdata/rawdata.imiss and rawdata/rawdata.lmiss. In the .imiss-file the N_MISS column denotes the number of missing SNPs, and the F_MISS column denotes the proportion of missing SNPs per individual.

The grey dashed line in Figure 8.1 indicates the mean call rate, while the red dashed line indicates the threshold we had determined above.

Figure 8.1: Per sample call rate.

8.3 Heterozygosity rate

To properly calculate heterozygosity rate and relatedness (identity-by-descent [IBD]) we need to do four things:

1. pre-clean the data to get a high-quality set,
2. of independent SNPs,
3. exclude long-range linkage disequilibrium (LD) blocks that bias with these calculations, and
4. exclude A/T and C/G SNPs as these may be ambivalent in interpretation when frequencies between cases and controls are close (MAF ± 0.45),
5. remove all non-autosomal SNPs.

You can find an up-to-date list of LD blocks you should exclude in these types of analyses here for the different genome builds. In this case we are using build 37. For the purpose of this book we included a file with these regions in the support-directory.

We will use the following settings:

• remove A/T and C/G SNPs with the flag --exclude dummy_project/all.atcg.variants.txt,
• call rate <1% with the flag --geno 0.10,
• Hardy-Weinberg Equilibrium (HWE) p-value > 1x10-3 with the flag --hwe 1e-3,
• and MAF>10% with the flag --maf 0.10,
• prune the data to only select independent SNPs (with low LD r^2) of one pair each with r^2 = 0.2 with the flags --indep-pairwise 100 10 0.2 and --extract rawdata/raw-GWA-data.prune.in,
• SNPs in long-range LD regions (for example: MHC chr 6 25.8-36Mb, chr 8 inversion 6-16Mb, chr17 40-45Mb, and a few more) with the flag --exclude range support/exclude_problematic_range.txt,
• remove non-autosomal SNPs with the flag --autosome.

First, get a list of A/T and C/G SNPs. Remember, the list of markers for this GWAS is noted in the .bim file. We can simply grep all the lines where the two alleles either have an A/T or C/G combination.

cat dummy_project/rawdata.bim | \
awk '($5 == "A" && $6 == "T") || ($5 == "T" && $6 == "A") || ($5 == "C" && $6 == "G") || ($5 == "G" && $6 == "C")' | awk '{ print $2, $1, $4, $3, $5, $6 }' \
> dummy_project/all.atcg.variants.txt

Second, clean the data and get a list of independent SNPs.

plink --bfile dummy_project/rawdata \
--autosome \
--maf 0.10 --geno 0.10 --hwe 1e-3 \
--indep-pairwise 100 10 0.2 \
--exclude range support/exclude_problematic_range.txt \
--make-bed --out dummy_project/rawdata.clean.temp

Please note, we have create a dataset without taking into account LD structure. Hence, the message ‘Pruned 0 variants from chromosome 1, leaving 19420.’ etc. In a dataset without any LD structure this flag --indep-pairwise 100 10 0.2 doesn’t actually work. However, with real-data you can use it to prune out unwanted SNPs in high LD.

Third, exclude the pruned SNPs. Note, how we include a file to exclude high-LD for the purpose of the practical.

plink --bfile dummy_project/rawdata.clean.temp \
--extract rawdata/raw-GWA-data.prune.in \
--make-bed --out dummy_project/rawdata.clean.ultraclean.temp

Fourth, remove the A/T and C/G SNPs.

plink --bfile dummy_project/rawdata.clean.ultraclean.temp \
--exclude dummy_project/all.atcg.variants.txt \
--make-bed --out dummy_project/rawdata.clean.ultraclean

Please note, this dataset doesn’t actually include this type of SNP, hence rawdata/all.atcg.variants.txt is empty! Again, you can use this command in real-data to exclude A/T and C/G SNPs.

Lastly, remove the temporary files.

rm -fv dummy_project/*.temp*

Finally, we can calculate the heterozygosity rate.

plink --bfile dummy_project/rawdata.clean.ultraclean --het --out dummy_project/rawdata.clean.ultraclean

This creates the file dummy_project/rawdata.clean.ultraclean.het, in which the third column denotes the observed number of homozygous genotypes, O(Hom), and the fifth column denotes the number of nonmissing genotypes, N(NM), per individual. We can now calculate the observed heterozygosity rate per individual using the formula (N(NM) - O(Hom))/N(NM).

Often there is a correlation between heterozygosity rate and missing data. Thus, we should plot the observed heterozygosity rate per individual on the x-axis and the proportion of missing SNP, that is the ‘SNP call rate’, per individuals on the y-axis (Figure 8.2).

Figure 8.2: Heterozygosity as a function of SNP call rate.

Examine the plot (Figure 8.2) to decide reasonable thresholds at which to exclude individuals based on elevated missing or extreme heterozygosity. We chose to exclude all individuals with a genotype failure rate >= 0.03 (vertical dashed line) and/or a heterozygosity rate ± 3 s.d. from the mean (horizontal dashed lines). Add the FID and IID of the samples failing this QC to the file named fail-imisshet-qc.txt.

How would you create this file?

If all is right, you’d have something like Table 8.2.

Table 8.2: Failed samples due to sample call rates and heterozygosity rate.

FID	IID
28	28
50	50
52	52
53	53
67	67
96	96
205	205
298	298
310	310
351	351
404	404
418	418
423	423
428	428
442	442
444	444
548	548
772	772
835	835
850	850
853	853
910	910
939	939
966	966
968	968
1,003	1,003
1,006	1,006
1,045	1,045
1,058	1,058
1,154	1,154
1,236	1,236
1,294	1,294
1,395	1,395
1,537	1,537
1,554	1,554
1,587	1,587
1,694	1,694
1,789	1,789
1,832	1,832
1,866	1,866
1,904	1,904
1,920	1,920
1,951	1,951
1,952	1,952
1,953	1,953
1,954	1,954
1,955	1,955
1,981	1,981
1,984	1,984
1,985	1,985
1,986	1,986

8.4 Relatedness

We calculate Identity-by-Descent (IBS) to identify duplicated and related samples. In Table 8.3 we show how much DNA is shared between individuals depending on their relation(Stapples J. 2014). IBS is measured by calculating pi-hat (\(\widehat{\pi}\)), which is in essence the proportion of the DNA that a pair of samples share. To calculate this, we needed this ultraclean dataset, without low-quality SNPs and without high-LD regions.

Table 8.3: Familial relations and % DNA shared.

Relatedness	%.DNA.sharing	IBD0	IBD1	IBD2
Monozygotic twins	±100%	0	1	0.00
Parents/child	±50%	0.25	0.5	0.25
Sibling	±50%	0.25	0.5	0.25
Fraternal twins	±50%	0.25	0.5	0.25
Grandparent/grandchild	±25%	0.5	0.5	0.00
Aunt/Uncle/Niece/Nephew	±25%	0.5	0.5	0.00
Half-sibling	±25%	0.5	0.5	0.00
First-cousin	±12.5%	0.75	0.25	0.00
Half first-cousin	±6.25%	0.875	0.125	0.00
First-cousin once removed	±6.25%	0.875	0.125	0.00
Second-cousin	±3.13%	0.9375	6.25E-2	0.00
Second-cousin once removed	±1.56%	0.96875	3.125E-2	0.00
Distantly related	<1.56%	varies	varies	0.00
Unrelated (includes relationships beyond the third degree)	<1.56%	1	0	0.00

PLINK calculates the inter-individual relatedness using the --genome function.

plink --bfile dummy_project/rawdata.clean.ultraclean --genome --out dummy_project/rawdata.clean.ultraclean

We can now identify all pairs of individuals with an IBD > 0.185. The code looks at the individual call rates stored in rawdata.imiss and outputs the IDs of the individual with the lowest call rate to ‘fail-IBD-QC.txt’ for subsequent removal (Table 8.4).

First, move to the dummy_project directory.

cd dummy_project

Now, execute this script - it should work just fine out-of-the-box.

perl ../scripts/run-IBD-QC.pl rawdata rawdata.clean.ultraclean

Go back one directory.

cd ..

Table 8.4: Failed IBD and callrate.

FID	IID
1,952	1,952
1,953	1,953
1,954	1,954
1,955	1,955
1,957	1,957
1,959	1,959
1,961	1,961
1,963	1,963
1,965	1,965
1,967	1,967
1,969	1,969
1,971	1,971
1,973	1,973
1,975	1,975

8.5 Ancestral background

Using a Principal Component Analysis (PCA) we can reduce the dimensions of the data, and project the “ancestral distances”. In other words, the principal component 1 (the first dimension) and principal component 2 (the second dimension) which will capture most of the variation in the data and represent how much each sample is alike the next. And when compared to a reference, you can deduce the ancestral background of each sample in your dataset. Of course this is relative: we will only know that a given sample is very much a like samples from a given population that exists today.

Nowadays we run such PCA against a large and diverse dataset containing many different populations. Old-school GWAS (pre-2009) would compare a dataset against HapMap 3, nowadays we prefer at a minimum the 1000G phase 3 populations. And in those ancient times the preferred software to run a PCA was Eigensoft which is a bit tricky to install (see Chapter 19), but nowadays PLINK provides the --pca-flag.

For the purpose of this practical primer we will run PCA using PLINK and its --pca-flag against an earlier version of 1000G, phase 1, which is slightly smaller and just as good to use.

8.5.1 1000G phase 1

We will project our data to the reference, in this example 1000G phase 1 (1000G), which includes individuals from 14 distinct global populations across 4 ‘super’-populations (Europeans [EUR], Africans [AFR], East-Asians [EAS], and Latin Americans [AMR]). In the real-world, using phase 1 may be just fine, but if you think your population evolved through extensive migration it’s probably best to use phase 3 data. In other words, the choice of reference is really depending on the dataset.

First, we will merge our data with 1000G. The alleles at each marker must be aligned to the same DNA strand to allow our data to merge correctly. Because not all SNPs are required for this analysis the A->T and C->G SNPs, which are more difficult to align, can be omitted.

8.5.1.1 Filter the 1000G data

First, we should get a list of relevant variants from our rawdata-dataset. We don’t need the other variants present in the 1000G dataset, right?

cat dummy_project/rawdata.bim | grep "rs" > dummy_project/all.variants.txt

Extract those from the 1000G phase 1 data.

plink --bfile ref_1kg_phase1_all/1kg_phase1_all --extract dummy_project/all.variants.txt --make-bed --out ref_1kg_phase1_all/1kg_phase1_raw

8.5.1.2 Filter A/T & C/G SNPs

As explained, the A/T and C/G SNPs are problematic, we want to exclude these too. So let’s get a list of A/T and C/G variants from 1000G to exclude - this may take a while.

cat ref_1kg_phase1_all/1kg_phase1_raw.bim | \
awk '($5 == "A" && $6 == "T") || ($5 == "T" && $6 == "A") || ($5 == "C" && $6 == "G") || ($5 == "G" && $6 == "C")' | awk '{ print $2, $1, $4, $3, $5, $6 }' \
> ref_1kg_phase1_all/all.1kg.atcg.variants.txt

Exclude those A/T and C/G variants in both datasets and at the same time filter to only retain high-quality data and exclude non-autosomal variants.

plink --bfile ref_1kg_phase1_all/1kg_phase1_raw --exclude ref_1kg_phase1_all/all.1kg.atcg.variants.txt -make-bed --out ref_1kg_phase1_all/1kg_phase1_raw_no_atcg

plink --bfile dummy_project/rawdata --exclude ref_1kg_phase1_all/all.1kg.atcg.variants.txt --make-bed --out dummy_project/rawdata_1kg_phase1_raw_no_atcg

8.5.1.3 Merging datasets

Try and merge the data.

plink --bfile dummy_project/rawdata_1kg_phase1_raw_no_atcg --bmerge ref_1kg_phase1_all/1kg_phase1_raw_no_atcg --make-bed --out dummy_project/rawdata.1kg_phase1

There probably is an error …

Error: 72 variants with 3+ alleles present.
* If you believe this is due to strand inconsistency, try --flip with
  dummy_project/rawdata.1kg_phase1-merge.missnp.
  (Warning: if the subsequent merge seems to work, strand errors involving SNPs
  with A/T or C/G alleles probably remain in your data.  If LD between nearby
  SNPs is high, --flip-scan should detect them.)
* If you are dealing with genuine multiallelic variants, we recommend exporting
  that subset of the data to VCF (via e.g. '--recode vcf'), merging with
  another tool/script, and then importing the result; PLINK is not yet suited
  to handling them.
See https://www.cog-genomics.org/plink/1.9/data#merge3 for more discussion.

So let’s flip some variants.

plink --bfile dummy_project/rawdata --exclude ref_1kg_phase1_all/all.1kg.atcg.variants.txt --flip dummy_project/rawdata.1kg_phase1-merge.missnp --make-bed --out dummy_project/rawdata_1kg_phase1_raw_no_atcg

Let’s try again and merge the data.

plink --bfile dummy_project/rawdata_1kg_phase1_raw_no_atcg --bmerge ref_1kg_phase1_all/1kg_phase1_raw_no_atcg --make-bed --out dummy_project/rawdata.1kg_phase1

There still is an error – there are a few multi-allelic variants present which PLINK can’t handle.

Error: 14 variants with 3+ alleles present.
* If you believe this is due to strand inconsistency, try --flip with
  dummy_project/rawdata.1kg_phase1-merge.missnp.
  (Warning: if the subsequent merge seems to work, strand errors involving SNPs
  with A/T or C/G alleles probably remain in your data.  If LD between nearby
  SNPs is high, --flip-scan should detect them.)
* If you are dealing with genuine multiallelic variants, we recommend exporting
  that subset of the data to VCF (via e.g. '--recode vcf'), merging with
  another tool/script, and then importing the result; PLINK is not yet suited
  to handling them.
See https://www.cog-genomics.org/plink/1.9/data#merge3 for more discussion.

Let’s just remove these multi-allelic variants.

plink --bfile dummy_project/rawdata_1kg_phase1_raw_no_atcg --exclude dummy_project/rawdata.1kg_phase1-merge.missnp --make-bed --out dummy_project/rawdata_1kg_phase1_raw_no_atcg_bi

After removing those pesky multi-allelic variants, we should be able to merge the data. We should take of the following:

• extract the pruned SNP-set (remember?), --extract rawdata/raw-GWA-data.prune.i,
• exclude non-autosomal variants, --autosome,
• and only keeping high-quality data, --maf 0.10 --geno 0.10 --hwe 1e-3,

plink --bfile dummy_project/rawdata_1kg_phase1_raw_no_atcg_bi \
--bmerge ref_1kg_phase1_all/1kg_phase1_raw_no_atcg \
--autosome \
--maf 0.10 --geno 0.10 --hwe 1e-3 \
--extract rawdata/raw-GWA-data.prune.in \
--make-bed --out dummy_project/rawdata.1kg_phase1.clean

Before we continue it’s best to clean up a bit of the mess.

rm -fv dummy_project/rawdata_1kg_phase1_raw_no_atcg_bi* ref_1kg_phase1_all/1kg_phase1_raw_no_atcg* dummy_project/rawdata.1kg_phase1.pruned* dummy_project/rawdata_1kg_phase1_raw_no_atcg* dummy_project/rawdata.1kg_phase1-merge* dummy_project/rawdata.1kg_phase1.bed dummy_project/rawdata.1kg_phase1.bim dummy_project/rawdata.1kg_phase1.fam dummy_project/rawdata.1kg_phase1.hh dummy_project/rawdata.1kg_phase1.log ref_1kg_phase1_all/1kg_phase1_raw.*

Now we have prepared our dataset with only high-quality SNPs that have few missing data, that are high-frequent, exclude problematic genomic ranges, and merged to the 1000G phase 1 reference dataset. Your output should look something like this:

PLINK v1.90b7.2 64-bit (11 Dec 2023)           www.cog-genomics.org/plink/1.9/
(C) 2005-2023 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to dummy_project/rawdata.1kg_phase1.clean.log.
Options in effect:
  --autosome
  --bfile dummy_project/rawdata_1kg_phase1_raw_no_atcg_bi
  --bmerge ref_1kg_phase1_all/1kg_phase1_raw_no_atcg
  --extract rawdata/raw-GWA-data.prune.in
  --geno 0.10
  --hwe 1e-3
  --maf 0.10
  --make-bed
  --out dummy_project/rawdata.1kg_phase1.clean

16384 MB RAM detected; reserving 8192 MB for main workspace.
2000 people loaded from dummy_project/rawdata_1kg_phase1_raw_no_atcg_bi.fam.
1092 people to be merged from ref_1kg_phase1_all/1kg_phase1_raw_no_atcg.fam.
Of these, 1092 are new, while 0 are present in the base dataset.
Warning: Multiple positions seen for variant 'rs3934834'.
Warning: Multiple positions seen for variant 'rs3737728'.
Warning: Multiple positions seen for variant 'rs6687776'.
Warning: Multiple chromosomes seen for variant 'rs1050301'.
Warning: Multiple chromosomes seen for variant 'rs4850'.
317476 markers loaded from dummy_project/rawdata_1kg_phase1_raw_no_atcg_bi.bim.
312239 markers to be merged from ref_1kg_phase1_all/1kg_phase1_raw_no_atcg.bim.
Of these, 14 are new, while 312225 are present in the base dataset.
312190 more multiple-position warnings: see log file.
Performing single-pass merge (3092 people, 308317 variants).
Merged fileset written to dummy_project/rawdata.1kg_phase1.clean-merge.bed +
dummy_project/rawdata.1kg_phase1.clean-merge.bim +
dummy_project/rawdata.1kg_phase1.clean-merge.fam .
308317 variants loaded from .bim file.
3092 people (1522 males, 1570 females) loaded from .fam.
3092 phenotype values loaded from .fam.
--extract: 49856 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 3078 founders and 14 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.994867.
299 variants removed due to missing genotype data (--geno).
--hwe: 13825 variants removed due to Hardy-Weinberg exact test.
2617 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
33115 variants and 3092 people pass filters and QC.
Phenotype data is quantitative.
--make-bed to dummy_project/rawdata.1kg_phase1.clean.bed +
dummy_project/rawdata.1kg_phase1.clean.bim +
dummy_project/rawdata.1kg_phase1.clean.fam ... done.

So in total there are 3,092 individuals, 1,522 males and 1,570 females, and 3,078 founders and 14 non-founders. The total genotyping rate is 99.5% and 33,115 variants are present.

8.5.2 Principal component analysis

Great, we’ve prepared our dummy project data and merged this with 1000G phase 1. Let’s execute the PCA using --pca in PLINK.

plink --bfile dummy_project/rawdata.1kg_phase1.clean --pca --out dummy_project/rawdata.1kg_phase1.clean

8.5.3 Plotting the PCA results

If all is peachy, you just succesfully ran PCA against 1000G phase 1. Using --pca we have calculated principal components (PCs), 20 in total by default, and we can now start plotting them. Let’s create a scatter diagram of the first two principal components, including all individuals in the file rawdata.1kg_phase1.clean.eigenvec (the first and second principal components are columns 3 and 4, respectively). We need to collect some per-sample information to color the points according to sample origin.

First we collect the results from the --pca, the dummy data phenotype information, and the reference population information.

Derive PC1 and PC2 thresholds so that only individuals who match the given ancestral population are included. For populations of European descent, this will be either the CEU or TSI 1000G individuals (Figure 8.3). Here, we chose to exclude all individuals with a first principal component score less than 0.0023.

Write the FID and IID of these individuals to a file called fail-ancestry-QC.txt.

cat dummy_project/rawdata.1kg_phase1.clean.eigenvec | \
awk '$3 < 0.0023' | awk '{ print $1, $2 }' > dummy_project/fail-ancestry-QC.txt

Choosing which thresholds to apply (and thus which individuals to remove) is not a straightforward process. The key is to remove those individuals with greatly divergent ancestry, as these samples introduce the most bias to the study. Identification of more fine-scale ancestry can be conducted by using less divergent reference samples (e.g., within Europe, stratification could be identified using the CEU, TSI (Italian), GBR (British), FIN (Finnish) and IBS (Iberian) samples from the 1,000 Genomes Project (http://www.1000genomes.org/)). Robust identification of fine-scale population structure often requires the construction of many (2–10) principal components.

Figure 8.3: PCA - Your data vs. 1000g.

8.6 Removing samples

Finally! We have a list of samples of poor quality or divergent ancestry, and duplicated or related samples. We should remove these. Let’s collect all IDs from our fail-*-files into a single file.

cat dummy_project/fail-* | sort -k1 | uniq > dummy_project/fail-qc-inds.txt

This new file should now contain a list of unique individuals failing the previous QC steps which we want to remove.

plink --bfile dummy_project/rawdata --remove dummy_project/fail-qc-inds.txt --make-bed --out dummy_project/clean_inds_data

Question: How many variants and samples are left? How many cases and how many controls did you loose?

8.7 The next step

Now that you filtered samples, we should turn our attention to step 2 of the QC for GWAS: identify SNPs of poor quality in Chapter @ref(gwas_basics_snp_qc).

9 Per-SNP QC

Now that we removed samples, we can focus on low-quality variants.

9.1 SNP call rates

We start by calculating the missing genotype rate for each SNP, in other words the per-SNP call rate.

plink --bfile dummy_project/clean_inds_data --missing --out dummy_project/clean_inds_data

Let’s visualize the results to identify a threshold for extreme genotype failure rate. We chose a callrate threshold of 3%, but it’s arbitrary and depending on the dataset, the data (visualization), and the number of samples (Figure 9.1).

Figure 9.1: Per SNP call rate.

9.2 Differential SNP call rates

There could also be differences in genotype call rates between cases and controls. It is very important to check for this because these differences could lead to spurious associations. We can test all markers for differences in call rate between cases and controls, or based on other criteria.

plink --bfile dummy_project/clean_inds_data --test-missing --out dummy_project/clean_inds_data

Let’s collect all the SNPs with a significantly different (P < 0.00001) missing data rate between cases and controls.

cat dummy_project/clean_inds_data.missing | awk '$5 < 0.00001' | awk '{ print $2 }' > dummy_project/fail-diffmiss-qc.txt

9.3 Allele frequencies

We should also get an idea on what the allele frequencies are in our dataset. Low frequent SNPs should probably be excluded, as these are uninformative when monomorphic (allele frequency = 0), or they may lead to spurious associations.

plink --bfile dummy_project/clean_inds_data --freq --out dummy_project/clean_inds_data

Let’s also plot these data. You can view the result below, and try it yourself (Figure 9.2).

Figure 9.2: Minor allele frequency.

9.3.1 A note on allele coding

Oh, one more thing about alleles.

PLINK codes alleles as follows:

A1 = minor allele, the least frequent allele A2 = major allele, the most frequent allele

And when you use PLINK the flag --freq or --maf is always relative to the A1-allele, as is the odds ratio (OR) or effect size (beta).

However, SNPTEST makes use of the so-called OXFORD-format, this codes alleles as follows:

A = the ‘other’ allele B = the ‘coded’ allele

When you use SNPTEST it will report the allele frequency as CAF, in other words the coded allele frequency, and the effect size (beta) is always relative to the B-allele. This means, CAF could be the MAF, or minor allele frequency, but this is not a given.

In other words, always make sure what the allele-coding of a given program, be it PLINK, SNPTEST, GCTA, et cetera, is! I cannot stress this enough. Ask yourself: ‘what is the allele frequency referring to?’, ‘the effect size is relative to…?’.

Right, let’s continue.

9.4 Hardy-Weinberg Equilibrium

Because we are performing a case-control genome-wide association study, we probably expect some differences in Hardy-Weinberg Equilibrium (HWE), but extreme deviations are probably indicative of genotyping errors.

plink --bfile dummy_project/clean_inds_data --hardy --out dummy_project/clean_inds_data

Let’s also plot these data. You can view the result below, and type over the code to do it yourself.

Figure 9.3: Hardy-Weinberg Equilibrium p-values per stratum.

9.5 Final SNP QC

We are ready to perform the final QC. After inspecting the graphs we will filter on a MAF < 0.01, call rate < 0.05, and HWE < 0.00001 (Figure 9.3), in addition those SNPs that failed the differential call rate test will be removed.

plink --bfile dummy_project/clean_inds_data --exclude dummy_project/fail-diffmiss-qc.txt --maf 0.01 --geno 0.05 --hwe 0.00001 --make-bed --out dummy_project/cleandata

9.6 A Milestone

Congratulations. You reached a very important milestone. Now that you filtered samples and SNPs, we can finally start the association analyses in Chapter @ref(gwas_testing).

10 Genome-Wide Association Study

Now that you have learned how to perform QC, you can easily run a GWAS and execute some downstream visualization and analyses. Let’s do this with a dummy dataset.

10.1 Exploring the data

Even though someone says that the QC was done, it is still wise and good practice to run some of the commands above to get a ‘feeling’ about the data. So let’s do this.

plink --bfile gwas/gwa --freq --out dummy_project/gwa

plink --bfile gwas/gwa --missing --out dummy_project/gwa

plink --bfile gwas/gwa --hardy --out dummy_project/gwa

Let’s visualize the results. First we should load in all the results.

Question: Load the data using R. [Hint: use and adapt the examples from the previous chapters.]

We can plot the per-stratum HWE p-values.

Question: Plot the per-stratum HWE p-values using R. [Hint: use and adapt the examples from the previous chapters.]

Figure 10.1: Per stratum HWE p-values.

We will want to see what the distribution of allele frequencies looks like.

Question: Plot the allele frequencies using R. [Hint: use and adapt the examples from the previous chapters.]

Figure 10.2: Minor allele frequencies.

We will want to identify samples that have poor call rates.

Question: Plot the per-sample call rates using R. [Hint: use and adapt the examples from the previous chapters.]

Figure 10.3: Per sample call rates.

We also need to know what the per SNP call rates are.

Question: Plot the per-SNP call rates using R. [Hint: use and adapt the examples from the previous chapters.]

Figure 10.4: Per SNP call rates.

10.2 Genetic models

A simple chi-square test of association can be done.

plink --bfile gwas/gwa --model --out gwas/data

Genotypic, dominant and recessive tests will not be conducted if any one of the cells in the table of case-control by genotype counts contains less than five observations. This is because the chi-square approximation may not be reliable when cell counts are small. For SNPs with MAFs < 5%, a sample of more than 2,000 cases and controls would be required to meet this threshold and more than 50,000 would be required for SNPs with MAF < 1%.

You can change this default behaviour by adding the flag --cell, e.g., we could lower the threshold to 3.

plink --bfile gwas/gwa --model --cell 3 --out gwas/data

Let’s review the contents of the results.

It contains 1,530,510 rows, one for each SNP, and each type of test (genotypic, trend, allelic, dominant, and recessive) and the following columns:

• chromosome [CHR],
• the SNP identifier [SNP],
• the minor allele [A1] (remember, PLINK always codes the A1-allele as the minor allele!),
• the major allele [A2],
• the test performed [TEST]:
  • GENO (genotypic association);
  • TREND (Cochran-Armitage trend);
  • ALLELIC (allelic as- sociation);
  • DOM (dominant model); and
  • REC (recessive model)],
• the cell frequency counts for cases [AFF],
• the cell frequency counts for controls [UNAFF],
• the chi-square test statistic [CHISQ],
• the degrees of freedom for the test [DF],
• and the asymptotic P value [P] of association.

Question: Do you know which model, i.e. TEST is most commonly used and reported? And why is that, do think?

10.3 Logistic regression

We can also perform a test of association using logistic regression. In this case we might want to correct for covariates/confounding factors, for example age, sex, ancestral background, i.e. principal components, and other study specific covariates (e.g. hospital of inclusion, genotyping centre etc.). In that case each of these P values is adjusted for the effect of the covariates.

When running a regression analysis, be it linear or logistic, PLINK assumes a multiplicative model. By default, when at least one male and one female is present, sex (male = 1, female = 0) is automatically added as a covariate on X chromosome SNPs, and nowhere else. The sex flag causes it to be added everywhere, while no-x-sex excludes it.

plink --bfile gwas/gwa --logistic sex --covar gwas/gwa.covar --out gwas/data

Let’s examine the results.

## [1] 918306      9

##      CHR       SNP      BP     A1   TEST NMISS     OR     STAT      P
##    <int>    <char>   <int> <char> <char> <int>  <num>    <num>  <num>
## 1:     1 rs3934834  995669      T    ADD  3818 1.0290  0.38120 0.7031
## 2:     1 rs3934834  995669      T    AGE  3818 1.0020  1.11800 0.2635
## 3:     1 rs3934834  995669      T    SEX  3818 1.0120  0.19090 0.8486
## 4:     1 rs3737728 1011278      A    ADD  3982 1.0190  0.38670 0.6990
## 5:     1 rs3737728 1011278      A    AGE  3982 1.0020  1.09800 0.2721
## 6:     1 rs3737728 1011278      A    SEX  3982 1.0060  0.09898 0.9212
## 7:     1 rs6687776 1020428      T    ADD  3915 0.9692 -0.33330 0.7389
## 8:     1 rs6687776 1020428      T    AGE  3915 1.0020  1.04000 0.2984
## 9:     1 rs6687776 1020428      T    SEX  3915 1.0150  0.23690 0.8127

Question: How come there are more lines in this file than there are variants?

If no model option is specified, the first row for each SNP corresponds to results for a multiplicative test of association. The C >= 0 subsequent rows for each SNP correspond to separate tests of significance for each of the C covariates included in the regression model. We can remove the covariate-specific lines from the main report by adding the hide-covar flag.

The columns in the association results are:

• the chromosome [CHR],
• the SNP identifier [SNP],
• the base-pair location [BP],
• the minor allele [A1],
• the test performed [TEST]: ADD (multiplicative model or genotypic model testing additivity),
  • GENO_2DF (genotypic model),
  • DOMDEV (genotypic model testing deviation from additivity),
  • DOM (dominant model), or
  • REC (recessive model)],
• the number of missing individuals included [NMISS],
• the OR relative to the A1, i.e. minor allele,
• the coefficient z-statistic [STAT], and
• the asymptotic P-value [P] of association.

We need to calculate the standard error and confidence interval from the z-statistic. We can modify the effect size (OR) to output the beta by adding the beta flag.

Question: Can you write down the mathematical relation between beta and OR?

10.4 Let’s get visual

Looking at numbers is important, but it won’t give you a perfect overview. We should turn to visualizing our results in Chapter 11.

11 GWAS visualisation

Data visualization is key, not only for presentation but also to inspect the results.

11.0.1 QQ plots

We should create quantile-quantile (QQ) plots to compare the observed association test statistics with their expected values under the null hypothesis of no association and so assess the number, magnitude and quality of true associations.

First, we will add the standard error, call rate, A2, and allele frequencies.

## Key: <SNP>
##           SNP   CHR        BP     A1 NMISS     OR    STAT       P     A2    MAF
##        <char> <int>     <int> <char> <int>  <num>   <num>   <num> <char>  <num>
## 1: rs10000010     4  21227772      C  3996 1.0420  0.9010 0.36760      T 0.4258
## 2: rs10000023     4  95952929      T  3957 0.9902 -0.2160 0.82900      G 0.4841
## 3: rs10000030     4 103593179      A  3991 0.9779 -0.3696 0.71170      G 0.1616
## 4:  rs1000007     2 237416793      C  4000 1.0180  0.3649 0.71520      T 0.3122
## 5: rs10000092     4  21504615      C  3963 0.9240 -1.6770 0.09354      T 0.3430
## 6: rs10000121     4 157793485      G  3919 0.9665 -0.7525 0.45170      A 0.4532
##    NCHROBS callrate
##      <int>    <num>
## 1:    7992  0.99900
## 2:    7914  0.98925
## 3:    7982  0.99775
## 4:    8000  1.00000
## 5:    7926  0.99075
## 6:    7838  0.97975

## [1] 306102     14

## # A tibble: 6 × 14
##   SNP        CHR     BP A1    A2      MAF callrate NMISS NCHROBS     BETA     SE
##   <chr>    <int>  <int> <chr> <chr> <dbl>    <dbl> <int>   <int>    <dbl>  <dbl>
## 1 rs10000…     4 2.12e7 C     T     0.426    0.999  3996    7992  0.0411  0.0457
## 2 rs10000…     4 9.60e7 T     G     0.484    0.989  3957    7914 -0.00985 0.0456
## 3 rs10000…     4 1.04e8 A     G     0.162    0.998  3991    7982 -0.0223  0.0605
## 4 rs10000…     2 2.37e8 C     T     0.312    1      4000    8000  0.0178  0.0489
## 5 rs10000…     4 2.15e7 C     T     0.343    0.991  3963    7926 -0.0790  0.0471
## 6 rs10000…     4 1.58e8 G     A     0.453    0.980  3919    7838 -0.0341  0.0453
## # ℹ 3 more variables: OR <dbl>, STAT <dbl>, P <dbl>

Let’s list the number of SNPs per chromosome. This gives a pretty good idea about the per-chromosome coverage. And it’s a sanity check: did the whole analysis run properly (we expect 22 chromosomes)?

Question: Why do the number of variants per chrosome (approximately) correlate with the chromosome number?

Question: Where are the data for chromosome X, Y and MT?

Let’s plot the QQ plot to diagnose our GWAS.

Figure 11.1: A QQ plot.

11.1 Manhattan plots

We also need to create a Manhattan plot to display the association test P-values as a function of chromosomal location and thus provide a visual summary of association test results that draw immediate attention to any regions of significance (Figure 11.2).

Figure 11.2: A manhattan plot.

11.2 Other plots

It is also informative to plot the density per chromosome. We can use the CMplot for that which you can find here. For now we just make these graphs ‘quick-n-dirty’, you can further prettify them, but you easily loose track of time, so maybe carry on.

Question: What do the grey spots on the density plot indicate?

This would lead to the following graphs.

Figure 11.3: SNP density of the association results.

Figure 11.4: A QQ plot including a 95% confidence interval (blue area) and genome-wide significant hits (red).

Figure 11.5: A regular manhattan plot. Colored by chromosome, suggestive hits are green, genome-wide hits are red. The bottom graph shows the per-chromosome SNP density.

Figure 11.6: A circular manhattan.

11.3 Interactive plots

You can also make an interactive version of the Manhattan - just because you can. The code below shows you how.

library(plotly)
library(dplyr)

# Prepare the dataset (as an example we use the data (gwasResults) from the `qqman`-package)
don <- gwasResults %>%

  # Compute chromosome size
  group_by(CHR) %>%
  summarise(chr_len=max(BP)) %>%

  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(chr_len)-chr_len) %>%
  select(-chr_len) %>%

  # Add this info to the initial dataset
  left_join(gwasResults, ., by=c("CHR"="CHR")) %>%

  # Add a cumulative position of each SNP
  arrange(CHR, BP) %>%
  mutate( BPcum=BP+tot) %>%

  # Add highlight and annotation information
  mutate( is_highlight=ifelse(SNP %in% snpsOfInterest, "yes", "no")) %>%

  # Filter SNP to make the plot lighter
  filter(-log10(P)>0.5)

# Prepare X axis
axisdf <- don %>% group_by(CHR) %>% summarize(center=( max(BPcum) + min(BPcum) ) / 2 )

# Prepare text description for each SNP:
don$text <- paste("SNP: ", don$SNP, "\nPosition: ", don$BP, "\nChromosome: ", don$CHR, "\nLOD score:", -log10(don$P) %>% round(2), "\nWhat else do you wanna know", sep="")

# Make the plot
p <- ggplot(don, aes(x=BPcum, y=-log10(P), text=text)) +

    # Show all points
    geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
    scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +

    # custom X axis:
    scale_x_continuous( label = axisdf$CHR, breaks= axisdf$center ) +
    scale_y_continuous(expand = c(0, 0)) +     # remove space between plot area and x axis

    # Add highlighted points
    geom_point(data=subset(don, is_highlight=="yes"), color="orange", size=2) +

    # Custom the theme:
    theme_bw() +
    theme(
      legend.position="none",
      panel.border = element_blank(),
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank()
    )
ggplotly(p, tooltip="text")

It will produce something like this.

Again, this is an example with dummy data - you can try to do it for our GWAS, but careful with the time. You can also choose to carry on.

You will encounter the above types of visualizations in any high-quality GWAS paper, because each is so critically informative. Usually, analysts of large-scale meta-analyses of GWAS will also stratify the QQ-plots based on the imputation quality (if your GWAS was imputed), call rate, and allele frequency (although that is rarely shared in publications, not even in supplemental material).

11.4 Stop playing around

Alright. It’s time to stop playing around and do a quick recap of what you’ve learned.

1. You learned how to convert datasets.
2. You learned how to execute sample QC and create diagnostic graphics
3. You learned how to do the same for SNP QC
4. You learned how to execute an association study given a dataset, covariates, and different assumptions regarding the genetic model.
5. You learned how to visualize results and played around with different visuals.

You should be ready for the real stuff. And if not, the next chapter will help you get ready: Chapter @ref(wtccc1_intro).

12 The Welcome Trust Case-Control Consortium

Now that you know your way around PLINK, bash and r and have done some basic quality control and association testing, you are ready for the real thing. We have prepared a real dataset: the first release of the Welcome Trust Case-Control Consortium (WTCCC) on coronary artery disease (CAD) and a control dataset used for that project.

12.1 Genotyping

The WTCCC1 data were genotyped using a chip from Affymetrix, nowadays part of ThermoFisher. As a brand Affymetrix still exists, but the chips aren’t made anymore. Unfortunately, most links to the old generation Affymetrix chips are borken, but you can still find some information about the 500K chip that was used for WTCCC1. It’s good practice to read up a bit on what chip was used, and what support materials are available.

12.2 The data

Before quality control the original data included:

• CAD cohort, n ± 2,000
• Healthy controls, from the UK 1958 birth control cohort, n ± 1,500 (we won’t use this)
• Healthy controls, from the UK National Blood Service, n ± 1,500.

12.3 Assignment

Your assignment in the next chapter (Chapter 13) is to do the following:

1. Explore the individual datasets by calculating some statistics and visualising these.
2. Merge the datasets in the folder wtccc1.
3. Calculate PCs using smartPCA.
4. Perform an association test using available covariates.
5. Visualize the results.
6. Identify independent SNPs.
7. Make regional association plots.

12.4 There you go

As I wrote, you are ready for the real stuff in Chapter 13.

13 WTCCC1: a GWAS on coronary artery disease (CAD)

As usual, we start by exploring the data in hand.

plink --bfile wtccc1/CADn1871_500Kb37fwd --bmerge wtccc1/UKBSn1397_500Kb37fwd --make-bed --out wtccc1/wtccc1 && \
plink --bfile wtccc1/wtccc1 --freq --out wtccc1/wtccc1 && \
plink --bfile wtccc1/wtccc1 --hardy --out wtccc1/wtccc1 && \
plink --bfile wtccc1/wtccc1 --missing --out wtccc1/wtccc1 && \
plink --bfile wtccc1/wtccc1 --test-missing --out wtccc1/wtccc1

cat wtccc1/wtccc1.missing | awk '$5 < 0.00001' | awk '{ print $2 }' > wtccc1/wtccc1-fail-diffmiss-qc.txt

Let’s investigate the HWE p-value in the whole cohort, and per stratum (cases and controls) with the code below.

This will result in Figure 13.1.

Figure 13.1: Stratified HWE p-values.

We should also inspect the allele frequencies. Note that by default PLINK (whether v0.7, v1.9, or v2.0) stores the alleles as minor (A1) and major (A2), and therefore --maf always calculates the frequency of the minor allele (A1).

This will result in Figure 13.2.

Figure 13.2: Minor allele frequencies.

There could be sample with very poor overall call rate, where for many SNPs there is no data. We will want to identify these samples and exclude them.

This will result in Figure 13.3.

Figure 13.3: Per sample call rate.

Lastly, we must inspect the per SNP call rate; we need to know if there are SNPs that have no data for many samples. We will want to identify such SNPs and exclude these.

This will result in Figure 13.4.

Figure 13.4: Per SNP call rate.

13.1 Quality control

Now that we have handle on the data, we can filter it.

Do you have any thoughts on that? Do you agree with the filters I set below? How would you do it differently and why?

plink --bfile wtccc1/wtccc1 --exclude wtccc1/wtccc1-fail-diffmiss-qc.txt --maf 0.01 --geno 0.05 --hwe 0.00001 --make-bed --out wtccc1/wtccc1_clean

13.2 Ancestral background

If these individuals are all from the United Kingdom, we are certain there will be admixture from other populations given UK’s history. Let’s project the WTCCC1 data on 1000G phase 1 populations.

We will face the same issue as before with our dummy dataset with respect to EIGENSOFT. So I created the data for you to skip to the Plotting PCA section immediately. Regardless, in the Preparing PCA and Running PCA sections I show you how to get there.

13.2.1 Preparing PCA

Filtering WTCCC1

For PCA we need to perform extreme clean.

plink --bfile wtccc1/wtccc1_clean --maf 0.1 --geno 0.1 --indep-pairwise 100 50 0.2 --exclude support/exclude_problematic_range.txt --make-bed --out wtccc1/wtccc1_temp

plink --bfile wtccc1/wtccc1_temp --exclude wtccc1/wtccc1_temp.prune.out --make-bed --out wtccc1/wtccc1_extrclean

rm -fv wtccc1/wtccc1_temp*

cat wtccc1/wtccc1_extrclean.bim | awk '{ print $2 }' > wtccc1/wtccc1_extrclean.variants.txt

cat wtccc1/wtccc1.bim | grep "rs" > wtccc1/all.variants.txt

Merging WTCCC1 with 1000G phase 1

Now we are ready to extract the WTCCC1 variants from the 1000G phase 1 reference

plink --bfile ref_1kg_phase1_all/1kg_phase1_all --extract wtccc1/all.variants.txt --make-bed --out ref_1kg_phase1_all/1kg_phase1_wtccc1

Extracting the A/T and C/G SNPs as well.

cat ref_1kg_phase1_all/1kg_phase1_wtccc1.bim | \
awk '($5 == "A" && $6 == "T") || ($5 == "T" && $6 == "A") || ($5 == "C" && $6 == "G") || ($5 == "G" && $6 == "C")' | awk '{ print $2, $1, $4, $3, $5, $6 }' \
> ref_1kg_phase1_all/all.1kg_wtccc1.atcg.variants.txt

plink --bfile ref_1kg_phase1_all/1kg_phase1_wtccc1 --exclude ref_1kg_phase1_all/all.1kg_wtccc1.atcg.variants.txt --make-bed --out ref_1kg_phase1_all/1kg_phase1_wtccc1_no_atcg

plink --bfile ref_1kg_phase1_all/1kg_phase1_wtccc1_no_atcg --extract wtccc1/wtccc1_extrclean.variants.txt --make-bed --out ref_1kg_phase1_all/1kg_phase1_raw_no_atcg_wtccc1

Finally we will merge the datasets.

plink --bfile wtccc1/wtccc1_extrclean --bmerge ref_1kg_phase1_all/1kg_phase1_raw_no_atcg_wtccc1 --maf 0.1 --geno 0.1 --exclude support/exclude_problematic_range.txt --make-bed --out wtccc1/wtccc1_extrclean_1kg

13.2.2 Running PCA

cp -v wtccc1/wtccc1_extrclean_1kg.bim wtccc1/wtccc1_extrclean_1kg.pedsnp
cp -v wtccc1/wtccc1_extrclean_1kg.fam wtccc1/wtccc1_extrclean_1kg.pedind

Now that the cleaning is done, we can execute the actual PCA.

perl ~/git/EIG/bin/smartpca.perl \
-i wtccc1/wtccc1_extrclean_1kg.bed \
-a wtccc1/wtccc1_extrclean_1kg.pedsnp \
-b wtccc1/wtccc1_extrclean_1kg.pedind \
-k 10 \
-o wtccc1/wtccc1_extrclean_1kg.pca \
-p wtccc1/wtccc1_extrclean_1kg.plot \
-e wtccc1/wtccc1_extrclean_1kg.eval \
-l wtccc1/wtccc1_extrclean_1kg.log \
-m 5 \
-t 10 \
-s 6.0 \
-w ref_1kg_phase1_all/1kg-pca-populations.txt

13.2.3 Plotting PCA

If all is peachy, you were able to run the PCA for the WTCCC1 data against 1000G phase 1. Using smartpca (you know, EIGENSOFT) we have calculated principal components (PCs) and we can now start plotting them. Let’s create a scatter diagram of the first two principal components, including all individuals in the file wtccc1_extrclean_1kg.pca.evec (the first and second principal components are columns 2 and 3, respectively). Use the data in column 4 to color the points according to sample origin.

Please note! You may have been able to make EIGENSOFT to work. So you may have to change “/ref_pca_wtccc1/wtccc1_extrclean_1kg.pca.evec” to “/wtccc1/wtccc1_extrclean_1kg.pca.evec” in the command below.

And we should visualize the PCA results: are these individuals really all from European (UK) ancestry?

We expect most individuals from the WTCCC to be 100% British, but a substantial group will have a different ancestral background as shown in the Figure 13.5 you just made.

Figure 13.5: PCA - WTCCC1 vs. 1000G

13.2.4 Removing samples

In a similar fashion as in the example gwas and rawdata datasets, you should consider to remove the samples below the threshold based on this PCA (Figure 13.5).

Go ahead, try that.

You’re code would be something like below:

cat wtccc1/wtccc1_extrclean_1kg.pca.evec | tail -n +2 | \
awk '$3 < 0.023' | awk '{ print $1 }' | awk -F":" '{ print $1, $2 }' > wtccc1/fail-ancestry-QC.txt

Next we filter these samples and get a final fully QC’d dataset.

plink --bfile wtccc1/wtccc1_clean --exclude wtccc1/fail-ancestry-QC.txt --make-bed --out wtccc1/wtccc1_qc

13.3 Association testing

Now that we have explored the data, we are ready for some simple association testing. However, it would be great to have some PCs to correct for. We can use PLINK for that too.

plink --bfile wtccc1/wtccc1_extrclean --exclude wtccc1/fail-ancestry-QC.txt --pca --out wtccc1/wtccc1_extrclean

Let’s add those PCs to the covariates-file.

echo "IID PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20" > wtccc1/wtccc1_qc.pca

cat wtccc1/wtccc1_extrclean.eigenvec | awk '{ print $2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$15,$16,$17,$18,$19,$20,$21,$22}' >> wtccc1/wtccc1_qc.pca

perl scripts/mergeTables.pl --file1 wtccc1/wtccc1_qc.pca --file2 wtccc1/wtccc1.covar --index IID --format NORM > wtccc1/wtccc1_qc.covar_pca

Exciting, now we are ready to perform a GWAS on CAD in WTCCC1.

plink --bfile wtccc1/wtccc1_qc --logistic sex --covar wtccc1/wtccc1_qc.covar_pca --out wtccc1/wtccc1_qc_log_covar_pca

After you ran the association analysis, you’re ready to process the data and take a first look at the results. First, we prepare the raw output.

You could visualize these results with the code below.

This would lead to the following graphs.

Figure 13.6: SNP density of the association results.

Figure 13.7: A QQ plot including a 95% confidence interval (blue area) and genome-wide significant hits (red).

Figure 13.8: A regular manhattan plot. Colored by chromosome, suggestive hits are green, genome-wide hits are red. The bottom graph shows the per-chromosome SNP density.

Figure 13.9: A circular manhattan.

13.4 Replication!

You reached an important milestone.

You recreated the work by the whole WTCCC1-team that took them years in just one afternoon!

Wow. Take a pause. And realize what you’ve done.

Back then there wasn’t much on analyses after a GWAS, nowadays there are many post-GWAS analyses methods. We will cover them in the next Chapter @ref(post_gwas).

14 Post-GWAS Analyses

A critical step in post-GWAS analysis is probably ‘mapping SNPs to genes’. It is critical, but it is also the most challenging. How do you even map SNPs to genes? What criteria to use? Should we take into account physical position? Or is it of interest that certain SNPs might influence downstream or upstream gene expression? And what of the fact that most loci discovered in GWAS are _inter_genic? What is the heritability of our trait? Are there any pleiotropic effects?

In the next few sections I deal with a couple of the downstream, post-GWAS analyses. We don’t have time to go practically into a few of these steps, but I do believe it is important to provide you with some pointers where to start.

14.1 Clumping

The Manhattan plot immediately draws your attention to the peaks above the genome-wide significance threshold. Clumping is the procedure in which we identify the independent hits, those top variants and the variants in linkage disequilibrium (their ‘LD buddies’) in the same genomic region (locus). The basic steps are as follows. First define a threshold above which to identify the top variant, usually this is the genome-wide significance threshold. Next, define the maximum p-value of association a variant may have with the trait of interest. Define the strength of the correlation that is allowed between the top variant and its LD buddies. And lastly define the size of region around the top variant to assess.

Clumping can easily be done in PLINK. And you can too with the dummy GWAS dataset in this additional chapter @ref(add_chapter_regional_plot). You don’t have to now, it’s just for fun.

14.2 Conditional analyses

The visualization of a GWAS is very appealing and draws attention to single peaks, and the most significant variant. And with clumping you can identify the independent signals, i.e. the top variants and their LD buddies in the locus. We implicitly assume that the top variant captures all the variation in the region by its linkage disequilibrium with an unknown causal variant, in other words it is ‘tagging’ the (un)measured potential causal variant in the region. Intuitively it is unlikely that a single causal variant is accounting for all the LD between the unknown causal variant and the measured (genotyped or imputed) variants at the locus. So, a naive focus on only the top variant ignores the fact that multiple causal variants exist and thus the variation captured by that single top variant underestimates the total variation that could be explained by that single region (locus).

Conditional analysis is a tool to identify the secondary variants associated with the trait of interest. This involves conditioning the association model (trait ~ Variant X + covariates) on the primary variant associated in that locus. A more comprehensive and general approach would be to condition on each variant in the whole genome starting with the top variants and stepwise selecting additional variants according to their conditional p-values. This strategy could potentially result in multiple top variants in the same locus. In other words, it could uncover different haplotypes that are causing the same association with the trait of interest.

This would be straightforward with super computers and full-access to individual level data. Unfortunately we often have neither.

Don’t despair, a method exists that fills this gap. With the program GCTA you can execute conditional analyses on GWAS summary statistics - the kind you just produced with the dummy data and with the WTCCC1 dataset - and a proper reference. You can read more on this in the paper.

For now, we will skip this part. I will add a whole chapter on this at some later stage.

14.3 Statistical finemapping

Statistical finemapping is closely related to conditional analysis. Where conditional analysis identifies secondary signals in a region associated to the trait of interest, statistical finemapping answers the question ‘which variants are likely causal to the trait?’. In more formal words, through statistical finemapping we identify the 95% credible set of causal variants. There are multiple methods and tools developed to get to this answer. I list a few and I encourage you to read up on these.

14.3.1 FINEMAP

This tool is very fast and versatile, and was developed by Christian Benner. It will identify the credible sets for each locus and calculates the posterior probabilities for each set. This is the more hands-on version of getting to credible sets.

14.3.2 LocusZoom

You could upload - privately or publicly - your data to my.locuszoom.org and obtain a list of causal variants with posterior probabilities. It uses a simple procedure to obtain the credible sets using these scripts. This is the more lazy version of getting a list of likely causal variants.

14.3.3 SuSIE

An laternate approach to credible set identification is through SuSIE. Through SuSIE you can identify credible sets under the assumption of multiple causal variants, whereas FINEMAP assumes a single causal variant.

For now, we will skip this part. I will add a whole chapter on this at some later stage.

14.4 FUMA: FUnctional Mapping and Annotation of GWAS

Researchers from the VUMC in Amsterdam have created an online tool that aids in mapping genes and function to GWAS: “Functional Mapping and Annotation of Genome-Wide Association Studies” a.k.a. FUMA. This online tool uses a variety of datasets and programs to prioritize genes and map these to associated loci.

We have covered some aspects of post-GWAS analyses, and a lot are covered by FUMA. Let’s try and annotate our WTCCC1 results. The assignment in this chapter 15 will be a bit more Do It Yourself.

14.5 Phenome-wide association study

A phenome-wide association study, or PheWAS, deals with assessing the association of top variants identified in GWAS with other traits. Through a PheWAS you’ll to get a notion whether on any of the top variants, the independent hits in your GWAS, have pleiotropic effects. In other words, whether your independent hits have effects on other traits too. This will paint a picture as it were about the role the genetic locus you identified plays in life: is it unique to your trait or does it affect other traits as well? What does it mean when your top variant near the gene FTO associates not only to BMI, but also to type 2 diabetes and coronary artery disease? PheWAS will not answer that question, but they do help in inventory all the other traits to which your hits are associated.

Again, we’ll get a bit hands-on and more Do It Yourself in the chapter that deals with this. Let’s jump over to Chapter 16.

14.6 Colocalization

From all of the above follows another intuitive question. Suppose a signal in your trait also shows association - to some extent - with other traits. Do these two signals than significantly overlap, more so than you would think based on chance? Colocalization deals with exactly this question.

We are not going to deal with colocalization, but a nice starting point is RACER about which you can read more in chapter @ref(add_chapter_regional_plot). I will add a whole chapter on this at some later stage.

14.7 Genetic correlation

Because of the underlying biology or because of the way they are measured or calculated, traits can be highly correlated - phenotypically. That is to say, when the levels of HDL are high, LDL is probably low, and so when you draw a correlation plot from some data obtained in a general population you’ll see a nice pattern. Biology may also cause traits to be genetically correlated: the same genetic variants influence multiple traits. In other words, if a variant is associated with higher levels of LDL-cholesterol, it may alternatively be associated with lower HDL-cholesterol, etc. Genetic correlation will aid in further understanding the relations between traits biologically: if there is a strong genetic correlation, it is may be due to the same biological pathways and so the traits are ‘linked’ through the same processes (maybe counter-acting processes).

Puzzle: what other phenomenom could cause a high but spurious genetic correlation?

We are not going to deal with genetic correlation. I will add a whole chapter on this at some later stage.

14.8 Causal inference: Mendelian Randomization

In observational studies we may find a strong association between a certain risk factor or biomarker and a disease. For instance, epidemiological studies show that high circulating cystatin C is associated with risk of cardiovascular disease (CVD), independent of creatinine-based renal function measurements(Laan S. W. 2016). However, residual confounding and reverse causality remain alternative explanations for the strong correlation between cystatin C and CVD, both of which are difficult to tease apart from traditional observational (epidemiological) studies.

Mendelian Randomization (MR) harnesses the properties of the genome to enable causal inference of a biomarker. Specifically, the invariant nature of the genome and the random distribution of alleles from parents to offspring at conception mean that genetic information is not influenced by disease status (reverse causality) and should be free from confounding by traditional risk factors(Laan S. W. 2016). So, genetic variation that modulates serum concentrations of cystatin C could serve as an instrumental variable to assess the effect of lifelong elevated concentrations of cystatin C on disease risk, independent of potential confounders(Laan S. W. 2016).

Finding the causal genes for complex disease is a primary objective of many researchers, because these genes are putative therapeutic targets. In this course we intend to find out whether Type 2 Diabetes (T2D) causes coronary artery disease and ischemic stroke.

There is an easy, quick-and-dirty way and a harder way.

14.8.1 MRBase

The quick-and-dirty way uses MRBase. It is based on the TwoSampleMR package and uses the MRInstruments package to load in all the genetic instruments. The creators of MRBase have curated hundreds of GWAS and molecular QTL studies, and prepared the data for easy use on the website and in the R packages.

14.8.2 TwoSampleMR

The hard way means getting your hands dirty and using the TwoSampleMR package yourself in R. You’ll learn all the steps you need to take into account to get to the answer of this question. You’ll create the diagnostic graphs and calculate the statistics.

14.8.3 Your choice

Your choice. In both cases you’ll learn how to do execute a MR, what to look out for in the statistics and diagnostic graphs, and how interpret the results. Your choice: Chapter @ref(mr_mrbase) or Chapter @ref(mr_twosamplemr).

In a future version I will aim to include a partial replication of the Cystatin C paper(Laan S. W. 2016)

14.9 Polygenic scores

Under the assumption of polygenicity, many variants with small effects contribute to the phenotypic variation in a trait or the risk to disease. Sample size and accurate phenotyping have a major impact on the power of a GWAS. It is estimated that 70% of variants that reach 10^-6 in an initial discovery GWAS will reach genome-wide significance with increasing sample sizes. The individual variants that do reach genome-wide significance the first or second time around do not explain all the phenotypic variation or residual disease risk.

To capture the polygenicity of a trait beyond individual genome-wide significant hits, we can calculate a polygenic score (PGS). There are several methods to calculate PGS, but essentially it comes down to multiplying the effect size at a given variant with the number of effect alleles a given individual carries. Depending on the method you can use the ‘hard-coded’ genotyped data, or the imputed data. In case of the latter the effect sizes are multiplied by the genotype probabilities or dosages, this ensures that the PGS takes into account the uncertainty intrinsic to imputed data.

Let’s consider the most intuitive PGS-method which is based on p-value thresholding, which was nicely applied in this classic paper. This method selects variants based on their p-value and than calculates the PGS. This can be done using increasingly liberal significance thresholds (p_T), for instance starting with all variants with p < 10^-8, next p < 10^-7, p < 10^-6, p < 10^-5, p < 10^-4, p < 10^-3, p < 0.05, p < 0.1, p < 0.2, p < 0.3, p < 0.4, p < 0.5. Under the assumption of polygenicity we expect an increasingly stronger correlation (r²) of the PGS at a certain p_T with the trait of interest. Another application of this method is provided by Van Setten(Setten J. 2015).

Lastly, there is a great explanatory website on how to calculate PGS, and it also shows one application of PGS (although much debate ensues about this particular application!).

14.10 What comes next

You come to the end of this practical primer. What is left are a summary (not written yet) of what you’ve learned and should take home. And some musings in the Epilogue (18 on this book and what the future holds.

15 Functional Mapping and Annotation of GWAS

15.1 Assignment

15.1.1 Tutorial

Go to the FUMA website, get an account, and study the online-tutorial.

15.1.2 Create the input

You will need to use the fwrite function in r to write the concatenated results of the WTCCC1 study (remember: wtccc1_assoc_compl).

Question: can you figure out the sample size of the WTCCC1 data you used?

Make sure you know what the column-names are in the file, you’ll need that for FUMA. You can use cat your_file | head to get the first 10 lines of the file.

You will need to compress the resulting output, this will make uploading to FUMA go faster. You can use the bash-program gzip for that. If you want to know what the options are for that program, you can use gzip -h.

15.1.3 Run FUMA - SNP to gene

Upload via the form on the FUMA-website. Since you’ve done the tutorial you are familiar with its options.

Select everything in the tabs Gene Mapping (positional mapping), Gene Mapping (eQTL mapping) but not GTEx v6 and GTEx v7, and Gene Mapping (3D Chromatin Interaction mapping) and leave the settings at Gene types and MHC region as-is. At MAGMA analysis set the MAGMA gene expression analysis to include all tissues, but GTEx v6 and GTEx v7.

Don’t forget to give your analysis a name.

This will take some time and so it’s a good moment to carry on with the rest of the practical or take a break, or study for the exam.

Questions

1. How many lead SNPs did we find?
2. What do the results of MAGMA (the gene-based test) look like and how many genes pass the threshold of multiple testing correction?
3. How many loci were mapped?
4. How many genes were physically located and how many were mapped to these loci?
5. Do you think all loci are ‘correct’, i.e. do you ‘believe’ all the signals looking at the mapping results? Why?
6. For what tissues are the signals enriched?
7. Are there any chromatin interactions discovered?

15.1.4 Run FUMA - Gene to function

Now that you mapped SNPs to genes, it’s time to go back to ‘My Jobs’. Select your job and perform GENE2FUNC.

Questions

1. What genes show the lowest expression across tissues?
2. And what genes the highest?
3. For what pathways are the signals enriched?
4. What molecular functions are mapped to the signals?

15.2 Some closing thoughts

FUMA is a great tool, but it comes with a caveat. It includes a couple of references of which it is not readily clear which variants are included - the authors do provide the codes used on Git, but still, you don’t know which variants precisely are filtered. That is key: perhaps the top variant you discovered is filtered in the reference. This means FUMA will not use it to map SNPs to genes, rather next best variant. This should be in high-LD - but, again, assumptions… And of course, the references used might not match your data well enough. So, my advice: use FUMA (why not be lazy rather than work hard?), but aware of such caveats as I described. All in all, I do think FUMA is very complete, intuitive, and it makes your work publication-ready because it creates just the right file-formats for you too (.png, .svg, .pdf, .jpeg).

That said, time to move on to inspect other phenotypes in relation to your findings in the next Chapter 16 or to return to the previous chapter on post-GWAS analyses @ref(post_gwas).

16 Phenome-Wide Association Study

Given the large biobanks available nowadays that have also genotyped the participants and collected a vast-array of information on them, it is possible to perform a Phenome-Wide Association Study.

There are several options in (Table 16.1).

Table 16.1: PheWAS Resources.

Site	Description
https://atlas.ctglab.nl	4,155 GWAS based on UK Biobank release 2
https://biobankengine.stanford.edu	UK Biobank based
http://pheweb.sph.umich.edu	University of Michigan EHR-based PheWAS
http://pheweb.sph.umich.edu/SAIGE-UKB/	1,403 GWAS UK Biobank based
http://phewas.mrbase.org	UK Biobank based

Assignment

1. Perform a PheWAS with a few of the resources and your favorite SNP from this tutorial.
2. Compare the different websites. What do you notice?
3. How is a PheWAS informative?

16.1 What’s next?

You now know to what other traits, risk factors and diseases our trait of interest (coronary artery disease) is correlated. An obvious next question would be: is this a causal correlation? We will cover this aspect in the next Chapter @ref(mendelian_randomization) or to return to the previous chapter on post-GWAS analyses @ref(post_gwas).

17 MRBase

You took the easy way out. No shame in that. Your choice. Right, so as I wrote, MRBase is very easy to handle. It is based on the TwoSampleMR package and uses the MRInstruments package to load in all the genetic instruments. The creators of MRBase have curated hundreds of GWAS and molecular QTL studies, and prepared the data for easy use on the website and in the R packages.

Using MRBase we intend to find out whether Type 2 Diabetes (T2D) causes coronary artery disease and ischemic stroke.

17.1 The end?

You are almost at the end. Time to return to the previous chapter on post-GWAS analyses @ref(post_gwas).

18 Epilogue

What started as a simple ‘let’s write a practical on how to do a GWAS’, escalated into this GitBook. My second book. My fifth child (as far as I know). I hope you found it to be useful and learned a bit.

That said, much can be improved and so don’t hesitate to contact me.

As with any proper epilogue I should not forget a heartfelt and honest thank you to the readers and users of this work. Gratitude also goes to my dear colleagues Charlotte Onland, Jessica van Setten, and Kristel van Eijk, and former colleague Sara Pulit who asked me back in 2017 to join the course as a lecturer. It has been a pleasure to work with and learn from you, and it has been a fun (and sometimes stressful) experience to teach this course.

19 Additional: EIGENSOFT

There used to be a time that the preferred software EIGENSOFT for Principal Component Analysis (PCA) was EIGENSOFT For many, it still is. However, EIGENSOFT is a bit challenging to make it work to say the least. You need to install some programs, and this is not always straightforward.

So, here’s the deal.

I will share the how-to for a macOS environment below in [EIGENSOFT] - this should work in a Linux environment too as macOS is UNIX-based. You can choose to try and make it work on your system (be it UNIX or macOS based) at home (or in the office).

However, I recommend that you use the --pca function which is present in PLINK v1.9 and up. This means you should probably simply skip this section and jump straight to Chapter 8.

19.1 Install homebrew

You need to install brew, the missing package-manager and accompanying packages that Apple didn’t provide.

/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

Next, check that everything is in order.

brew doctor

19.2 Install missing packages

Right, now that you’ve done that, you’re ready to install gsl and openblas.

brew install gsl
brew install openblas

You may also require llvm for EIGENSOFT to work.

brew install llvm

19.3 Installing EIGENSOFT

I am still sharing the code you’ll need - you could try this on your system.

mkdir -v $HOME/git
cd $HOME/git
git clone https://github.com/DReichLab/EIG.git
cd EIG/src
make
make install

20 Answers to the questions

Throughout the book you encountered questions related to some steps and output that you generated. Here you can review the answers for each of these questions - of course, for some different answers are possible.

[ANSWERS NEED TO BE WORKED OUT]

20.1 Chapter 3.1: GWAS Basics

Question: Can you think off what the ‘11,440 heterozygous haploid genotypes present’ represent?

Question: Can you think of other scenarios in which you may want to extend the check on differential missingness beyond a check between cases and controls? * For instance ‘genotyping platform’, or ‘hospital of inclusion’, if you think this might influence the genotyping experiment technically.

20.2 Chapter 3.4: GWAS Basics - Association testing

Question: Load the data using R. [Hint: use and adapt the examples from the previous chapters.]

library("data.table")

COURSE_loc = "~/Desktop/practical" # getwd()

gwas_HWE <- data.table::fread(paste0(COURSE_loc, "/dummy_project/gwa.hwe"))
gwas_FRQ <- data.table::fread(paste0(COURSE_loc, "/dummy_project/gwa.frq"))
gwas_IMISS <- data.table::fread(paste0(COURSE_loc, "/dummy_project/gwa.imiss"))
gwas_LMISS <- data.table::fread(paste0(COURSE_loc, "/dummy_project/gwa.lmiss"))

Question: Plot the per-stratum HWE p-values using R. [Hint: use and adapt the examples from the previous chapters.]

library("ggpubr")

gwas_HWE$logP <- -log10(gwas_HWE$P)

ggpubr::gghistogram(gwas_HWE, x = "logP",
                    add = "mean",
                    add.params = list(color = "#595A5C", linetype = "dashed", size = 1),
                    rug = TRUE,
                    # color = "#1290D9", fill = "#1290D9",
                    color = "TEST", fill = "TEST",
                    palette = "lancet",
                    facet.by = "TEST",
                    bins = 50,
                    xlab = "HWE -log10(P)") +
  ggplot2::geom_vline(xintercept = 5, linetype = "dashed",
                      color = "#E55738", size = 1)
ggplot2::ggsave(paste0(COURSE_loc, "/dummy_project/gwas-hwe.png"),
       plot = last_plot())

Question: Plot the allele frequencies using R. [Hint: use and adapt the examples from the previous chapters.]

ggpubr::gghistogram(gwas_FRQ, x = "MAF",
                    add = "mean", add.params = list(color = "#595A5C", linetype = "dashed", size = 1),
                    rug = TRUE,
                    color = "#1290D9", fill = "#1290D9",
                    xlab = "minor allele frequency") +
  ggplot2::geom_vline(xintercept = 0.05, linetype = "dashed",
                      color = "#E55738", size = 1)
ggplot2::ggsave(paste0(COURSE_loc, "/dummy_project/gwas-freq.png"),
       plot = last_plot())

Question: Plot the per-sample call rates using R. [Hint: use and adapt the examples from the previous chapters.]

gwas_IMISS$callrate <- 1 - gwas_IMISS$F_MISS

ggpubr::gghistogram(gwas_IMISS, x = "callrate",
                    add = "mean", add.params = list(color = "#595A5C", linetype = "dashed", size = 1),
                    rug = TRUE, bins = 50,
                    color = "#1290D9", fill = "#1290D9",
                    xlab = "per sample call rate") +
  ggplot2::geom_vline(xintercept = 0.95, linetype = "dashed",
                      color = "#E55738", size = 1)
ggplot2::ggsave(paste0(COURSE_loc, "/dummy_project/gwas-sample-call-rate.png"),
       plot = last_plot())

Question: Plot the per-SNP call rates using R. [Hint: use and adapt the examples from the previous chapters.]

gwas_LMISS$callrate <- 1 - gwas_LMISS$F_MISS

ggpubr::gghistogram(gwas_LMISS, x = "callrate",
                    add = "mean", add.params = list(color = "#595A5C", linetype = "dashed", size = 1),
                    rug = TRUE, bins = 50,
                    color = "#1290D9", fill = "#1290D9",
                    xlab = "per SNP call rate") +
  ggplot2::geom_vline(xintercept = 0.95, linetype = "dashed",
                      color = "#E55738", size = 1)
ggplot2::ggsave(paste0(COURSE_loc, "/dummy_project/gwas-snp-call-rate.png"),
       plot = last_plot())

Question: Do you know which model, i.e. TEST is most commonly used and reported? And why is that, do think?

Question: How come there are more lines in this file than there are variants?

Question: Can you write down the mathematical relation between beta and OR?

20.3 Chapter 3.5: GWAS visualization

Question: Why do the number of variants per chrosome (approximately) correlate with the chromosome number?

Question: Where are the data for chromosome X, Y and MT?

Question: What do the grey spots on the density plot indicate?

20.4 Chapter 5.4: FUMA

Question: can you figure out the sample size of the WTCCC1 data you used?

Questions

1. How many lead SNPs did we find?
2. What do the results of MAGMA (the gene-based test) look like and how many genes pass the threshold of multiple testing correction?
3. How many loci were mapped?
4. How many genes were physically located and how many were mapped to these loci?
5. Do you think all loci are ‘correct’, i.e. do you ‘believe’ all the signals looking at the mapping results? Why?
6. For what tissues are the signals enriched?
7. Are there any chromatin interactions discovered?

Questions

1. What genes show the lowest expression across tissues?
2. And what genes the highest?
3. For what pathways are the signals enriched?
4. What molecular functions are mapped to the signals?

21 Licenses and disclaimers

21.1 Copyright

This book and all its material (“content”) is protected by copyright under Dutch Copyright laws and is the property of the author or the party credited as the provider of the content. You may not copy, reproduce, distribute, publish, display, perform, modify, create derivative works, transmit, or in any way exploit any such content, nor may you distribute any part of this content over any network, including a local area network, sell or offer it for sale, or use such content to construct any kind of database. You may not alter or remove any copyright or other notice from copies of the content on this website. Copying or storing any content except as provided above is expressly prohibited without prior written permission of the author or the copyright holder identified in the individual content’s copyright notice. For permission to use this content, please contact the author.

21.2 Disclaimer

The content contained herein is provided only for educational and informational purposes or as required by Dutch law. The author attempted to ensure that content is accurate and obtained from reliable sources, but does not represent it to be error-free. The author may add, amend or repeal any text, procedure or regulation, and failure to timely post such changes to this book shall not be construed as a waiver of enforcement. The author does not warrant that any functions on this website or the contents and references herein will be uninterrupted, that defects will be corrected, or that this website or the contents and references will be free from viruses or other harmful components. Any links to third party information on the author’s website are provided as a courtesy and do not constitute an endorsement of those materials or the third party providing them.

21.3 Images and data used

I took the at-most care to use refer to the original works and data sources where needed. Likewise, all the images c.q. figures are either produced specifically for this book, or I took them from Unsplash to brighten up the book. If you feel I made a mistake and your work should be properly referenced, please don’t hesitate to contact me.

These are the images from Unsplash listed here in no particular order.

• papers_on_wall - https://unsplash.com/photos/open-book-lot-Oaqk7qqNh_c
• women_behind_macbook - https://unsplash.com/photos/woman-using-macbook-pro-with-person-in-white-top-bPVM4nOy0Rg
• woman_working_on_code - https://unsplash.com/photos/woman-wearing-black-t-shirt-holding-white-computer-keyboard-YK0HPwWDJ1I
• licenses - https://unsplash.com/photos/book-lot-on-black-wooden-shelf-zeH-ljawHtg

21.4 Copyright

Copyright 1979-2024. All rights reserved. Sander W. van der Laan | s.w.vanderlaan [at] gmail.com | https://vanderlaanand.science. Published with bookdown.

22 Colophon

The 2022 and 2024 editions of this book were produce in RStudio and with the bookdown package. Below a listing of installed programs and libraries, the operating system, and their specific versions.

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.3.3 (2024-02-29)
##  os       macOS Sonoma 14.5
##  system   x86_64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/New_York
##  date     2024-04-04
##  pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package           * version date (UTC) lib source
##  askpass             1.2.0   2023-09-03 [2] CRAN (R 4.3.0)
##  bookdown          * 0.38.1  2024-03-26 [2] Github (rstudio/bookdown@50a1c1e)
##  bslib               0.6.2   2024-03-22 [2] CRAN (R 4.3.2)
##  cachem              1.0.8   2023-05-01 [2] CRAN (R 4.3.0)
##  chromote            0.2.0   2024-02-12 [1] CRAN (R 4.3.2)
##  cli                 3.6.2   2023-12-11 [2] CRAN (R 4.3.0)
##  colorspace          2.1-0   2023-01-23 [2] CRAN (R 4.3.0)
##  crayon              1.5.2   2022-09-29 [2] CRAN (R 4.3.0)
##  crul                1.4.0   2023-05-17 [2] CRAN (R 4.3.0)
##  curl                5.2.1   2024-03-01 [2] CRAN (R 4.3.2)
##  data.table          1.15.4  2024-03-30 [1] CRAN (R 4.3.2)
##  digest              0.6.35  2024-03-11 [2] CRAN (R 4.3.2)
##  evaluate            0.23    2023-11-01 [2] CRAN (R 4.3.0)
##  fastmap             1.1.1   2023-02-24 [2] CRAN (R 4.3.0)
##  flextable         * 0.9.5   2024-03-06 [1] CRAN (R 4.3.2)
##  fontBitstreamVera   0.1.1   2017-02-01 [2] CRAN (R 4.3.0)
##  fontLiberation      0.1.0   2016-10-15 [2] CRAN (R 4.3.0)
##  fontquiver          0.2.1   2017-02-01 [2] CRAN (R 4.3.0)
##  formatR           * 1.14    2023-01-17 [2] CRAN (R 4.3.0)
##  gdtools             0.3.7   2024-03-05 [2] CRAN (R 4.3.2)
##  gfonts              0.2.0   2023-01-08 [2] CRAN (R 4.3.0)
##  glue                1.7.0   2024-01-09 [2] CRAN (R 4.3.0)
##  htmltools           0.5.8   2024-03-25 [2] CRAN (R 4.3.2)
##  httpcode            0.3.0   2020-04-10 [2] CRAN (R 4.3.0)
##  httpuv              1.6.15  2024-03-26 [2] CRAN (R 4.3.2)
##  jquerylib           0.1.4   2021-04-26 [2] CRAN (R 4.3.0)
##  jsonlite            1.8.8   2023-12-04 [2] CRAN (R 4.3.0)
##  kableExtra        * 1.4.0   2024-01-24 [1] CRAN (R 4.3.2)
##  knitr             * 1.45    2023-10-30 [1] CRAN (R 4.3.0)
##  later               1.3.2   2023-12-06 [2] CRAN (R 4.3.0)
##  lifecycle           1.0.4   2023-11-07 [2] CRAN (R 4.3.0)
##  magrittr            2.0.3   2022-03-30 [2] CRAN (R 4.3.0)
##  mime                0.12    2021-09-28 [2] CRAN (R 4.3.0)
##  munsell             0.5.1   2024-04-01 [1] CRAN (R 4.3.2)
##  officer             0.6.5   2024-02-24 [2] CRAN (R 4.3.2)
##  openssl             2.1.1   2023-09-25 [2] CRAN (R 4.3.0)
##  processx            3.8.4   2024-03-16 [2] CRAN (R 4.3.2)
##  promises            1.2.1   2023-08-10 [2] CRAN (R 4.3.0)
##  ps                  1.7.6   2024-01-18 [2] CRAN (R 4.3.0)
##  R6                  2.5.1   2021-08-19 [2] CRAN (R 4.3.0)
##  ragg                1.3.0   2024-03-13 [2] CRAN (R 4.3.2)
##  Rcpp                1.0.12  2024-01-09 [2] CRAN (R 4.3.0)
##  rlang               1.1.3   2024-01-10 [2] CRAN (R 4.3.0)
##  rmarkdown         * 2.26.1  2024-03-26 [2] Github (rstudio/rmarkdown@ee69d59)
##  rstudioapi          0.16.0  2024-03-24 [2] CRAN (R 4.3.2)
##  sass                0.4.9   2024-03-15 [2] CRAN (R 4.3.2)
##  scales              1.3.0   2023-11-28 [2] CRAN (R 4.3.0)
##  sessioninfo         1.2.2   2021-12-06 [2] CRAN (R 4.3.0)
##  shiny               1.8.1   2024-03-26 [2] CRAN (R 4.3.2)
##  stringi             1.8.3   2023-12-11 [2] CRAN (R 4.3.0)
##  stringr             1.5.1   2023-11-14 [2] CRAN (R 4.3.0)
##  svglite             2.1.3   2023-12-08 [1] CRAN (R 4.3.0)
##  systemfonts         1.0.6   2024-03-07 [2] CRAN (R 4.3.2)
##  textshaping         0.3.7   2023-10-09 [2] CRAN (R 4.3.0)
##  tinytex           * 0.50    2024-03-16 [2] CRAN (R 4.3.2)
##  uuid                1.2-0   2024-01-14 [2] CRAN (R 4.3.0)
##  viridisLite         0.4.2   2023-05-02 [2] CRAN (R 4.3.0)
##  webshot           * 0.5.5   2023-06-26 [1] CRAN (R 4.3.0)
##  webshot2          * 0.1.1   2023-08-11 [1] CRAN (R 4.3.0)
##  websocket           1.4.1   2021-08-18 [1] CRAN (R 4.3.0)
##  xfun                0.43    2024-03-25 [2] CRAN (R 4.3.2)
##  xml2                1.3.6   2023-12-04 [2] CRAN (R 4.3.0)
##  xtable              1.8-4   2019-04-21 [2] CRAN (R 4.3.0)
##  yaml                2.3.8   2023-12-11 [2] CRAN (R 4.3.0)
##  zip                 2.3.1   2024-01-27 [2] CRAN (R 4.3.2)
## 
##  [1] /Users/slaan3/Library/R/x86_64/4.3/library
##  [2] /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────

References

Anderson C. A., et al. 2010. “Data QC in Genetic Case-Control Association Studies.” https://www.ncbi.nlm.nih.gov/pubmed/21085122.

Laan S. W., et al. van der. 2016. “Cystatin c and Cardiovascular Disease, a Mendelian Randomization Study.” https://doi.org/10.1016/j.jacc.2016.05.092.

Laurie C. C., et al. 2010. “Quality Control and Quality Assurance in Genotypic Data for Genome-Wide Association Studies.” https://www.ncbi.nlm.nih.gov/pubmed/20718045.

Purcell S., et al. 2007. “PLINK, a Tool Set for Whole-Genome Association and Population-Based Linkage Analyses.” https://www.ncbi.nlm.nih.gov/pubmed/17701901.

Setten J., et al. van. 2015. “Serum Lipid Levels, Body Mass Index, and Their Role in Coronary Artery Calcification: A Polygenic Analysis.” https://doi.org/10.1161/CIRCGENETICS.114.000496.

Stapples J., et al. 2014. “PRIMUS: Rapid Reconstruction of Pedigrees from Genome-Wide Estimates of Identity by Descent.” https://doi.org/10.1016/j.ajhg.2014.10.005.